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Q Abstract 



The issue of computing (co)homology generators of a cell complex is gaining a piv- 
otal role in various branches of science. While this issue can be rigorously solved in 
polynomial time, it is still overly demanding for large scale problems. Drawing inspi- 
ration from low- frequency electrodynamics, this paper presents a physics inspired 
algorithm for first cohomology group computations on three-dimensional complexes. 
The algorithm is general and exhibits orders of magnitude speed up with respect to 
competing ones, allowing to handle problems not addressable before. In particular, 
when generators are employed in the physical modeling of magneto-quasistatic prob- 
lems, this algorithm solves one of the most long-lasting problems in low- frequency 
computational electromagnetics. In this case, the effectiveness of the algorithm and 
its ease of implementation may be even improved by introducing the novel concept 
of lazy cohomology generators. 
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1 Introduction 



The availability of unprecedented computing power and efficient numerical 
methods produced a dramatic increase in applications of computational (co)homology [1] 
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(i.e. the computation of generators of the (co)homology group.) (Co) homology, 
in fact, has been already shown to be essential in unexpected areas of sci- 
ence, ranging from computer aided design (CAD) for feature detection [2J, 
parametrization and mesh generation [3] , shape analysis and pattern recogni- 
tion [4J, to sensors networks [5] and robot motion planning [6], medicine and 
biology [7J and quantum chemistry [iSj. Focusing on physics, (co)homology gen- 
erators have been used for example to detect the chaotic behavior in sampled 
experimental data [9], [ID], [EE], in quantum information theory [12], string 
theory p3] and electromagnetism [H], [15], [33], [TF], [18], [32]. 

The problem addressed in this paper is the computation of the first cohomol- 
ogy group generators for three-dimensional combinatorial manifolds, being the 
computation of the zeroth and second (co)homology groups already satisfac- 
torily solved in literature. 

Cohomology generators over integers — unlike the real and complex ones — can 
be rigorously computed in polynomial time by finding the Smith normal form 
(SNF) pQ of the coboundary matrix. However, this approach is computation- 
ally not attractive because the best implementation of the SNF exhibits a 
hyper-cubical complexity [20] • Sparse matrix data structures and the reduc- 
tion of the input complex [13], [17] before the SNF computation is run have 
been used to speed up the computation to a point that can be used for prac- 
tical problems. Recently, a novel algorithm called Thinned Current Technique 
(TCT) has been introduced in [21]. This algorithm exhibits orders of magni- 
tude speed up on typical problems with respect to other algorithms presented 
in literature and it is easy to implement. Nonetheless, its main limitation is 
that when the complex is not skeletonizable — i.e. cannot be homotopically re- 
tracted to a graph — this algorithm uses again the technique presented in [T7J . 
For a comprehensive survey on other algorithms proposed in literature refer 
to pE]. 

To introduce a physics inspired algorithm for first cohomology group compu- 
tations, let us focus on the interplay between (co)homology, discrete Hodge 
decomposition [22], [E] arid physical modeling that appears when consider- 
ing problems whose definition of potentials is not straightforward. One of the 
most studied examples where this happens occurs in low-frequency electrody- 
namics. Electromagnetic phenomena are governed by Maxwell's laws [23] and 
material constitutive relations. For slowly time- varying fields, whose change in 
magnetic field energy is dominant and electromagnetic wave propagation can 
be ignored, it is typical to brake the symmetry of Maxwell's laws by neglecting 
the displacement current in the Ampere-Maxwell's equation [23]. Using this 
magneto-quasistatic (MQS) approximation, the magnetic field is irrotational 
in the insulators thanks to Ampere's law. The fact that the insulating region is 
in most cases not simply connected prevent it to be exact. The consequence is 
that a magnetic scalar potential cannot be introduced naively in the insulating 
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regions. Yet, using a scalar potential is tempting since formulations based on 
it are computationally much more efficient than the ones using the classical 
magnetic vector potential. 

How to define a magnetic scalar potential in non simply connected domains 
has drawn a considerable effort in the computational electromagnetics commu- 
nity in the last twenty-five years. A connection of this issue with (co)homology 
theory has been advocated many years ago [21], [2]. Despite most scientists 
keep using other heuristic and sometimes patently incorrect definition of po- 
tentials and related algorithms (see references and counter-examples in [19]). 
it has been recently shown that the first cohomology group generators over 
integers of the cell complex modeling the insulating region are needed to make 
the problem well defined [IE], [TS]. Introductory material on this subject com- 
prising an informal introduction on algebraic topology, how to model physical 
variables as cochains with complex coefficients and how to rephrase Maxwell's 
laws in algebraic form, can be found in [16], [19], [21 J . 

The quest for an algorithm for first cohomology group computation that is 
both general and exhibits a linear average complexity is still open. This is 
surprising since the research on this issue has been pushed forward by many 
leading software houses having at least part of the core business in solving 
MQS problems as MagSoft Corp., Ansoft Corp.-ANSYS Inc., Vector Fields 
L.t.d., CST Ag., and Comsol Inc.. The fact that this issue has been considered 
unsolved for so many years indicates that computing cohomology generators 
quickly (and correctly) is not straightforward. Also the implementation com- 
plexity may affect negatively the technology transfer. Developing a simple and 
fast algorithm would enable to embed it in the next-generation of electromag- 
netic Computer-Aided Engineering (CAE) softwares. 

This article fills this gap by exploiting a novel physics inspired approach 
to compute cohomology generators suitable for physical modeling called the 
Dlotko-Specogna (DS) algorithm. Moreover, the novel concept of lazy coho- 
mology generators is introduced to speed up and simplify the implementation 
when cohomology generators are employed in computational physics. 

The paper is structured as follows. Section[2]is a mild introduction to (co)homology 
theory that may be skipped by readers familiar with this topic. In Section [3] 
the role of (co)homology theory in low-frequency electrodynamics is recalled. 
The Dlotko-Specogna (DS) algorithm and the novel concept of lazy cohomology 
generators are introduced in Section |4j In Section [5] some numerical experi- 
ments are presented to compare the novel algorithm to other state-of-the-art 
algorithms in terms of efficiency and robustness. Finally, in Section [6j the 
conclusions are drawn. 
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Fig. 1. On the left (upper and lower), representants of homology generators of the 
annulus. On the right (upper and lower), representants of cohomology generators of 
the annulus. 

2 Mild introduction to algebraic topology 

In this section a mild introduction to algebraic topology is provided. For 
a rigorous one please consult [25J. Let us consider a discretization /C of a 
given three-dimensional space as cell complex (more precisely, a regular CW- 
complex |25j.) For simplicity, one may think about a simplicial complex. K. is 
a combinatorial manifold if the link of every vertex is a sphere or disk. The 
link of a vertex v £ K, consists of all elements s G JC that do not contain v 
that are faces of higher dimensional elements containing v. 

Homology and its dual cohomology theory are mathematical tools to describe 
"holes" of a given space in a rigorous way. 

As an example, let us consider the two-dimensional simplicial complex K of the 
annulus represented in the Fig. [T] The zero-dimensional holes are defined as the 
connected components of /C. In the example, /C is formed by a single connected 
component. Clearly, there is one (one-dimensional) hole in the annulus. This 
hole can be surrounded with a one-dimensional oriented cycle, as in Fig. [T] 
upper left. This kind of cycle represents a generator of the first homology 
group of the annulus. 



The concept of holes may be generalized to higher dimensions. In the example, 
there are no holes of dimensions two and higher. Yet, two-dimensional holes 
for three-dimensional complexes are voids as the ones that may be encircled by 
a ball contained in K, (i.e. cavities of the complex.) In this case, the collection 
of oriented two-dimensional cells on the surface of the ball represents a second 
homology generator of the complex. The number of i-dimensional holes is often 
referred to as the zth Betti number /3»(/C). 

In homology theory some cycles are considered the same. If two i-dimensional 
cycles (with orientation) form a boundary of a set made of (i + 1) -dimensional 
elements, then we say that they are in the same homology (equivalence) class, 
or simply that they are equal. Two cycles in the same homology class are 
represented in Fig. [T] lower left. If some cycle alone is a boundary, then we say 
that it is homologically trivial. 

More formally, a chain c of a complex K, is a formal sum of elements of K, with 
coefficients (in this paper we consider integer or complex coefficients.) Chain c 
is a cycle if its boundary dc vanishes, c is homologically nontrivial if it is not a 
boundary. The group of z-dimensional cycles is denoted by Z^JC). The group 
of i-dimensional boundaries is denoted by Bi()C). The ith homology group is 
the quotient Hi{K) = Z^K) / B^K). By Z) and C) we denote the 

integer and complex homology groups, respectively. 

Homology has a dual theory called cohomology. In case of complexes embedded 
in a tree-dimensional space, a straightforward duality between homology and 
cohomology generators exists. Let us fix a zth homology generator c. The 
dual cohomology generator c* is a set of ith dimensional elements having the 
following property: every cycle in the class of c needs to cross c*. Moreover, 
once the discrete analogous of integration of c* on c is considered by means 
of the dot product (c*,c) between the vectors representing the coefficients of 
c* and c, the result is always 1 (for more details consult [26].) A possible 
cohomology generator in the considered example can be visualized in Fig. [T] 
upper right. Analogously as in the case of homology, in cohomology two c* and 
c* are considered equal if there exist a set s of (i — l)-dimensional elements 
such that c* and c*' are common coboundary of elements of s, see Fig. [I] 
bottom right for an example. 

More formally, a cochain with integer coefficients is a map from the group of 
chains to integers. A cochain c* is a cocycle if its coboundary 5c* vanishes. 5 
is the coboundary operator [25J defined with the Generalized Stokes Theorem 
{5c*, c) = (c*,dc). c* is cohomologically nontrivial if it is not a coboundary. 
The group of i-dimensional cocycles of the complex K is denoted by Z % {K\ 
whereas the group of z-dimensional coboundaries is denoted as B l ()C). The ith 
cohomology group is the quotient H 1 (1C) = Z % {1C) / B % (K). By _£P(/C,Z) and 
H l (fC,C) we denote integer and complex cohomology groups, respectively. By 
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a (co)homology basis we mean a set of (co)homology generators that span the 
(co)homology group. 

One may also consider the so-called relative (co) homology group of a complex 
JC modulo a sub-complex /Co C JC in which /Q is forgotten. The relative 
(co)homology groups are denoted by i?j(/C, /Co) and fP(/C,/Co). The details 
can be found in [25] or, more informally, in [26], [18J. 

The dual complex JC pQ, [HI Section 3] is obtained from /C by using the 
bary centric subdivision. Let us define the dual cell complex /C = D(JC) in the 
following way: 

(1) For every polyhedron t G JC, n — D(t) is defined as the bary center of £. 

(2) For every triangle / G JC that is a common face of polyhedra ti,t 2 G /C, 
e = D(f) is defined as the sum of a segment of line joining the barycenter 
of / with D(t\) and a segment of line joining the barycenter of / with 
D(t 2 ). 

(3) For every edge e G JC let /i, . . . , f n G JC be the triangles incidental to e. 
f = D(e) is then defined as a (set-theoretical) sum \J™ =1 conv[_B(e), Z)(/j)], 
where -B(e) denotes the barycenter of the edge e and conv[-] the convex 
hull. 

(4) For every node n G JC let ei, . . . , e n G JC be the edges incidental to n. 
i = D{n) is the volume bounded by -D(ei), . . . , D(e n ). 



3 (Co) homology and low-frequency electrodynamics 

Let us cover the topologically trivial computational domain by a conformal 
polyhedral mesh JC which is a regular CW-complex [25]. Two sub-complexes 
JC C and JC a of JC are introduced that contain mesh elements belonging to 
the conducting and insulating regions, respectively. Let us define the poten- 
tials in JC a for a harmonic analysis of a MQS problem formulated by using a 
magnetic scalar potential, as the T-f2 formulation [TB], |19j . [2"T] . The alge- 
braic Ampere's law in the insulating region is enforced on every 2-cell with 
SF = I = 0, where I is the complex-valued electric current 2-cochain and F is 
the magneto-motive force (m.m.f.) complex- valued 1-cochain. Thanks to the 
discrete Hodge decomposition, the 1-cocycle F G Z 1 (/C a ,C) can be expressed 
as a linear combination of a basis of the 1st cohomology group H l (JC a ,C) 
plus a 1-coboundary _B 1 (/C a ,C). The 1-coboundary _B 1 (/C a ,C) is provided by 
taking the O-coboundary of a complex-valued 0-cochain f2 whose coefficients 
represent the magnetic scalar potential sampled on mesh nodes. Since JC a is 
embedded in M 3 , the (co)homology groups are torsion free [lj and the basis of 
i/ 1 (/C a , C) can be obtained from a basis of H 1 ()C a , Z) where the elements of Z 
are treated as elements of C [25]. Then, the nonlocal (i.e. applied not locally 
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Fig. 2. On the left, the independent current flowing in a torus. On the right, the 
1-cycle Cj £ Ci(/C ) and a 2-chain Sj such that Cj = dsj. 



on each 2-cell as SF = I, but on an arbitrary 2-chain) algebraic Ampere's 
law [IB] . [T9] (F,Cj) = (I, Sj) (see Fig. [2]) implicitly holds for any 1-cycle 
Cj G Ci(/C Q ), with Cj = dsj, by setting 

F = sn+ ]T (i) 



where (., .) denotes the dot product, {h- J }^i 1 are the representatives of the 
1st cohomology group iJ 1 (/C a ,Z) generators and ^i(JC a ) is the 1st Betti num- 
ber of K-a- When it is not confusing, by cohomology generators we refer to both 
the cohomology classes and their representatives. Physically, the can 
be interpreted as a set of independent currents [16j, [IS] flowing in the branches 
of the conductors /C c . Fig. [2] shows the independent current flowing in a solid 
2-dimensional torus. For a n-fold 2-dimensional solid torus, there are n inde- 
pendent currents. 

Let us fix the generators {h J }^ 1 Ca ' ) of ff 1 (/C a ,Z). There exists a set of cy- 
cles {ci}i=i being a Z) basis such that (h J , q) = 5y holds [2S], [2S]- 
Since /C is homologically trivial, there exist 2-chains {sj}£I\ a G C^/C) whose 
boundaries are the {cj}^ . Their restrictions {<jj}fl^ Ca ' ) to K, c form a H 2 (IC C , d}C c , Z) 
basis [21], see Fig. |2j on the right. The independent currents are exactly 
the currents linked in nonlocal algebraic Ampere's law [16J by the 1-cycles 
{ Cl }tf a) in such a way that the j-th independent current can be defined with 
(F, Cj) = ij = (I, Sj) = (I, <jj). Then, the currents linked by any other cycle in 
K, a can be obtained by linear combinations of the independent currents. 
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Fig. 3. On the left, the two cohomology generators for a 2-torus intersecting in the 
thicker edge. On the right, the support of the thinned current t 1 corresponding to 
c 1 . 

4 The Dlotko-Specogna (DS) algorithm 

Usually, cohomology generators are found first [16] , [19] and the current distri- 
bution inside /C c is a result of the MQS simulation. In principle, taking inspi- 
ration from physics, one may do it the other way around: The j-th 1-cocycle 
h ? can be computed by imposing a unity current in Oj and a zero current in 

Pi ( fC ^ 

the other sigmas, where the {crj}^ are H2(fC c ,dfC c ,Z) generators. If a set 
of currents {t 1 , . . . , t^ 1 ^")} that fulfill these constraints is constructed, finding 
the cohomology basis is just a matter of solving /?i(/C a ) linear systems 5 la? = P 
and setting a zero coefficient of the output cochains in JC c \K, a . Surprisingly, 
these systems can be solved in most cases by back-substitutions only, without 
using any integer linear solver and even a sparse matrix data structure, with 
the Extended Spanning Tree Technique (ESTT) [27], [21]. Intuitivelly ESTT 
algorithm is a way to extend a cohomology generator of a space 5* C X to a 
co cycle in X. 

The core of this physics-based approach to cohomology computation is how 
to obtain such currents. The key observation is that the 1-cocyles {h J } :/ I 1 
do not depend on how the unity currents are distributed inside K c . Moreover, 
if each unity current flows in a ring whose sections are "single 2-cells", then 
the constraint can be trivially imposed by assigning a unity current (the sign 
depends on incidence) to those 2-cells. The resulting 2-cocycles in /C c are 
called thinned currents [21J. An example of thinned current in a solid 2-torus 
is provided on the right of Fig. [3] In TCT algorithm [21] . a physically-based 
method to construct thinned currents is employed. A thinning is applied on /C c 
in such a way that the conductors become a "single 3-cell thick." The thinned 
conductors can be viewed on the dual complex as a graph representing the 
skeleton of K c . By computing independent cycles on that graph, by orienting 
them and by considering the 2-cells that are dual to the dual 1-cells in the 
graph, the thinned currents are found. This approach presents two problems: 
on the one hand, this approach does not work for conductors that do not 
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homotopically retract to a graph, on the other hand finding the skeleton is 
time consuming since all elements of /C c have to be processed. Here is when 
the DS algorithm, summarized as follows, comes is: 

(1) Find the n combinatorial 2-manifolds that represent the connected com- 
ponents Ci, . . . ,C n of dK, c (the external boundary dJC is excluded from 
/C c .) This part requires 0(card(K,)) time. 

(2) Compute the 1st cohomology B x {Ci, Z) generators c 1 , . . . , c 2g for each i- 
th connected component of <9/C c , where g is the genus of C;. This can be 
done in linear time worst-case complexity 0(card(dK. c ) g) with the graph- 
theoretic algorithm presented in [28] (see App.[A|for a short presentation.) 
For an example, see on the left of Fig. [3} 

(3) For every connected component Ci, find the thinned currents t 1 , . . . ,t 2g 
corresponding to c 1 , ... , c 2g in 0(card(dJC c ) g) with the following algo- 
rithm: 

for each 1-cell E with nonzero coefficient in c' 

for each 2-cell T G /C c with E in the boundary 

(t i ,T)+ = c E K{T,E) ] 
The value of the cochain t on a cell E is (t, E), whereas k(A, B) denotes 
the incidence between cells A and B, see [25]. Initially, set (t*,T) = 
for all 2-cells T G IC C . For an example, see on the right of Fig. [3j The 
demonstration of correctness of this approach is presented in App. [B} 

(4) For every connected component Ci, solve the integer systems ^h- 7 = t J , 
j G {1, . . . , 2g}, to find the 2g 1-cocycles h 1 , . . . , h 2g in K a . This can 
be performed without solving any system by a simultaneous application 
of the ESTT algorithm [27], [21]. Simultaneous means that the ESTT al- 
gorithm is applied to all t 1 , . . . ,t 2g thinned currents at the same time. 
Algorithmically this can be easily achieved by changing a real number to 
a vector of 2g real numbers in the ESTT algorithm. 

(5) For every connected component Ci, store the restrictions of h 1 , . . . , h 2g to 
K a . The average computational effort required is 0(card{fC) g). 

If the genus g is bounded by a constant 0(1), as it happens always in practical 
problems, the average complexity of the DS algorithm is linear 0(card(K.)). It 
is also purely graph-theoretic, so straightforward to implement. 

It is easy to prove that the 1-cocycles obtained by the DS algorithm span 
if 1 (/C a , Z). However, the obtained cocycles are not a basis of H 1 (JC a , Z), since 
the number of obtained generators is twice the cardinality of its own basis. 
In other words, some cocycle is a linear combination of the others. For this 
reason we refer to these redundant cocycles produced by the DS algorithm as 
lazy cohomology generators. 

To obtain a if 1 (/C a ,Z) cohomology basis with the DS algorithm, only the 
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cohomology generators on d)C c that are extended with the ESTT algorithm to 
coboundaries in JC a have to be used. They may be found by a change of basis 
obtained by computing linking numbers and a SNF of a matrix (see App. O 
for details.) Due to conceptual and implementation difficulties, one welcomes 
the possibility to bypass this additional step. 



Another contribution of this paper is to show that in computational physics 
avoid partitioning the i/ 1 (Ci,Z) generators is not only possible but even more 
convenient. This is realized by employing directly the lazy cohomology genera- 
tors in the physical modeling, something that has been never considered in the 
literature. Lazy cohomology generators are employed in a MQS formulation, 
for example the T-f2 [T5], [TB], [T9], as if they were a set of standard if 1 (/C a , Z) 
generators. Namely, a nonlocal Faraday's equation [16], [19] is written on the 
support of the j-th cohomology generator as (U, dhj) = — i uj (<&, hj), where U 
is the electro-motive force 1-cochain on the dual complex, <fr is the magnetic 
flux 2-cochain on the dual complex and hj = Dih^), D being the dual map [1] 
that maps elements of the original complex to elements of the dual complex. 
Lazy cohomology generators contain a H 1 (/C a Z) basis and generators that are 
dependent to the basis. Therefore, adding the dependent equations is not a 
problem considering that the system of equations to solve is already overde- 
termined. This is due to the fact that algebraic Faraday's equations [IB] . [19] 
enforced in K c are also dependent. Even though a full-rank system may be ob- 
tained by a tree-cotree gauging (i.e. set the electric vector potential on a tree 
of 1-cells in /C c to zero), it is widely known that with iterative linear solvers it 
is much more efficient to use an ungauged [29] formulation. What is important 
from the modeling point of view is that even if the potentials are not unique, 
the fields are. Therefore, the use of linearly dependent cocycles in the physical 
modeling does not introduce either any inconsistency in the formulation of the 
boundary value problem or any penalties in the computational time employed 
by the simulation due, for example, to a hypothetical increase of the condition 
number of the linear system matrix or to the use of twice as many cohomology 
generators as needed. 



As a final observation about lazy generators, let us now consider a lazy gen- 
erator belonging to the trivial class of if 1 (/C a ,Z). Given an arbitrary 1-cycle 
c G Zx(/C a ), the dot product of the lazy generator with c is zero. Therefore, 
trivial generators verify trivially the nonlocal algebraic Ampere's law and, in 
this case, the current ij does not represent the current linked by the dual ho- 
mology generator. Therefore, the value of the independent current relative to 
a trivial generator is not unique and it is determined by the solution of the 
system of equations. This is not surprising, since the independent current in 
this case does not have a physical meaning. 
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Table 1 



Time required (in seconds) for cohomology computation with various algorithms. 
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5 Numerical experiments 



The DS algorithm has been integrated into the research software CDICE [3T?] 
implemented in Fortran 90. Three competing algorithms presented in [TB], [T7j 
and [2T] are considered. Two computers are used to the computations: Intel 
Core 2 Duo T7700 2.4GHz laptop with 4GB of RAM and 64 GB of RAM 
and Intel Xenon E7-8830 2.13 GHz processors computer. As the mesh size 
increases, standard cohomology computations end up in failures due to having 
exceeded the memory limit of the laptop. When the laptop runs out of memory, 
the large computer is used. The DS algorithm has been executed on the laptop 
up to five millions of tetrahedra without encountering any problem, which 
shows that it is quite economical in terms of memory usage. 

Table [l] reports the time required (in seconds) for cohomology computation 
by the various algorithms. The timings obtained on the large computer are 
given in brackets. The Table presents also the time that DS algorithm requires 
for computing a standard cohomology basis (see App. [C] for implementation 
details.) 

The algorithm is going to be exploited to solve MQS problems arising in fusion 
engineering and design, in engineering and optimization of electromagnetic 
devices and the analysis of features of magnetic fields generated by current- 
carrying thick knots, see Fig. |4} 

We remark that the DS algorithm may be used also in applications outside com- 
putational physics where only the complex K a is available. In that case, in fact, 
one may employ an efficient mesh generator as TetGen (http : //www . tetgen . org 
to produce the mesh (without adding Steiner points |31j) of the complement 
with respect to a box containing K a . 
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Fig. 4. On the left, a complicated thick knot and the first cohomology group gener- 
ator of its complement. For clarity, the dual faces in the support of hj = D{hP) are 
shown. On the right, the current density obtained by the MQS solver. 

6 Conclusions 

The DS algorithm, even though is general and straightforward to implement, 
outperforms all competing state-of-the-art algorithms for first cohomology 
group computations. The time required for computing cohomology genera- 
tors with the DS algorithm is so limited that it allows to remove what has 
been considered for more than twenty-five years the main simulation bottle- 
neck for low-frequency electrodynamics problems. Therefore, we expect the 
DS algorithm to be embedded in the next-generation electromagnetic CAE 
softwares. 



A Cohomology generators of 2-manifolds 

In this section the algorithm to compute H 1 (d)C c ) generators presented in [28] 
is recalled. For simplicity, we consider 2-manifolds without boundary only. The 
extension to the general case can be found in [28J . 

By the primal skeleton of dlC c we mean the graph consisting of all the ver- 
tices and edges in dlC c . By the dual skeleton of dJC c we mean a graph whose 
vertices are the 2-cells in dlC c and an edge is put between two vertices iff. the 
corresponding faces in dfC c share an edge in dlC c . We want to point out that 
edges of both the primal and dual skeletons correspond to edges of dfC c . 

Let us fix a spanning tree T of the primal skeleton. Let us also fix a spanning 
tree T' of dual skeleton. We assume, that T and T' do not share edges of dfC c . 
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In [28] and [32] it is shown that the number of edges in <9/C c that are neither in 
T not in T' is the first Betti number of <9/C c . Moreover, the Hi(dJC c ) generators 
are the cycles closed in T by those edges, whereas the // 1 (9/C C ) generators are 
the cycles closed in T' by those edges. 



The idea of the procedure is presented in Figure A.l 





Fig. A.l. Upper left, the standard triangulation of a torus. Opposite sides are identi- 
fied. Upper right, tree on primal (solid bold), and dual (dotted bold) skeleton. With 
double bold, the two edges not belonging to one of the trees are depicted. Lower left, 
a cohomology generator closed by the first edge. Lower right, cohomology generator 
closed by the second edge. 



In order to obtain the coefficients of the cocycle, a simple procedure that 
orients the cycle is used. Let v and w be the vertices of the edge that close 
the cycle. With a BFS strategy [33] a distance function on the tree (or dual 
tree) from v is built as long as w is not reached. Then, a path in a tree from w 
to v is found by following the decreasing values of the defined function. The 
obvious details are left for the reader. 
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B Thinned currents construction 



In this section we show that the output t of the algorithm to obtain thinned 
currents from if 1 (9/C c ) is indeed a cochain (this is a necessary and sufficient 
condition for the ESTT algorithm termination [27].) To do this, first we need 
to remind one of Massey's equation |25j which is a property of any regu- 
lar CW-complex. Let T be a 3-dimensional cell having 1-dimensional cell E 
in boundary Then there exist exactly two 2-dimensional cells T x and T 2 in 
boundary of T both having E in their boundary. Moreover, the incidence 
indices satisfy the following equation (Theorem IX. 7. 2 in [25J): 



In order to show that £ is a cocycle, we have to show that St = 0. This is 
equivalent to showing that for every 3-dimensional cell T, (St,T) = 0. Of 
course, due to the the algorithm to obtain thinned currents from H l (dlC c ), 
this expression may be nonzero only for 3-dimensional cells in /C c that have 
some edge(s) from in input 1-cocycle c in boundary. To show that in this case 
(St,T) = 0, we need to consider the following cases: 

(1) T has only one edge E in the boundary such that (c, E) ^ 0. Then the 
two 2-cells T\ and T 2 from the Massey's equation will be nonzero in the 
cochain t (we have (t, Ti) = (c, E)k(Ti, E) and (t, T 2 ) = (c, E)k(T 2 , E) re- 
spectively.) From the Massey's equation, we know that k(T, Ti)k(Ti,E) + 
k(T, T 2 )k(T 2 , E) = 0. After multiplying this equation by (c, E) we get 



that prove this case. 
(2) T has two edges E\ and E 2 in the boundary such that (c,E\) ^ and 
(c, E 2 ) 7^ 0. We assume that E\ and E 2 are two constitutive edges in 
cocycle c, therefore they need to share a 2-cell T 3 in <9/C c . Since c is a 
cocycle, we have 



k(K, T x )k(T x , E) + k(K, T 2 )k(T 2 , E) = 0. 



k(T, T\)k{T\, E) (c, E) + k{T, T 2 )k{T 2 , E) (c, E) = 0. 



Therefore, 



= k(T,T 1 ){k(T 1 ,E)(c,E)) + 
K {T,T 2 ){k{T 2 ,E){c,E)) = k(T, Ti)(t, Ti) + 
K(T,T 2 )(t,T 2 ) = (St,T) 



(Et, c)k(T 3 , £i) + (E 2 , c)k(T 3 , E 2 ) = 0. 



(B.l) 



Moreover, we consider the following Massey's equations 



k(T, TO^Ti, E,) + k(T, T 3 )k(T 3 , E ± ) = 0, 



(B.2) 
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k(T, T 2 )n(T 2 , E 2 ) + k(T, T 3 )k{T 3 , E 2 ) = 0. 



(B.3) 



From (B.l), we have 



o = -«(t,t 3 )« j e; 1 ,c)«(T3, j e; 1 ) + (e 2 , c )k(t 3 ,e 2 )) = 

-k(T, T 3 )k(T 3 , E 1 ))(E 1 , c) + (-k(T, T 3 )k{T 3 , E 2 ))(E 2 , c). 



From (B.2) and (B.3), we have 

(-«(T, T 3 )k(T 3 , c) + (-«(T, T 3 )/<T 3 , £ 2 )) (£ 2 , c) = 

k(T, Ti)(k(Ti, E 1 )(E 1 ,c)) + k(T, T 2 )(k(T 2 , E 2 )(E 2 , c» = 
«(T,T 1 )(t,r 1 ) + K(T,T 2 )(t,r 2 > = <5t,T>, 

that prove this case. 
(3) If T has more than two edges in the boundary that have nonzero value in 
cochain c, then the conclusion follows from the inductive application of 
the argument for two edges. The obvious details are left for the reader. 



C Finding cocycles in dlC c that extend to a cohomology basis of JC a 



A technique to partition the homology generators on a combinatorial 2-manifold 
into two classes — the ones that bound in /C c , and the ones that bound in K a — 
is proposed in [31]. This algorithm is based on the computation of linking 
numbers between all possible pairs of generators and subsequent SNF com- 
putation of a small integer matrix containing the computed linking numbers. 
Linking numbers are defined for disjoint 1-cycles only, so a pre-processing has 
been applied in [34J on surface generators to perturb them in such a way no 
intersections are present between any pair of generators. For this purpose, 
the submerged cycles and shifted surface cycles have been found. This is not 
completely trivial to do algorithmically and cycles have to be "straightened" 
before [33]. The complexity of the linking number computation is quadratic 
with the number of edges in dK. c and a computationally costly interval arith- 
metic [35] package has to be used to rigorously compute linking numbers 
without the risk of errors due to the finite precision of real numbers [35] , [36] , 
especially on coarse meshes. 

In this paper, to find the change of basis to find cohomology generators in K. a 
we use the same idea presented in [34J, but finding linking numbers between 
paths on the dual complex pQ, [16]. This approach, contrarily to [33], does 
not require any extra computations since the surface cycles are obtained by 
considering the dual 1-cells that are dual to 1-cells in the support of the 
cohomology generators, while the submerged cycles are defined simply as the 
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dual 1-cells that are dual to 2-cells in the support of thinned currents. We 
want to remark that for this purpose the 2-cells belonging to the support of 
thinned currents have to be in order. They can be easily ordered in linear time 
and the technical details are left for the reader. 

About practical implementation of the change of basis, the algorithm pre- 
sented in [36] is used to compute linking numbers and the interval arithmetic 
library for Fortran 90 presented in [37] is employed for interval arithmetic 
computations. 
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